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Abstract. Wc investigate the Josephson effect for a setup with two lattice quantum 
wires featuring Majorana zero energy boundary modes at the tunnel junction. In the 
weak-coupling, the exact solution reproduces the perturbative result for the energy 
containing a contribution ~ ±cos(0/2) relative to the tunneling of paired Majorana 
fcrmions. As the tunnel amplitude g grows relative to the hopping amplitude w, the 
gap between the energy levels gradually diminishes until it closes completely at the 
critical value g,. = V^w. At this point the Josephson energies have the principal values 
= 2(7 v^i'^ COS [0/6 + 27r(m — l)/3], where m = —1,0,1 and a = ±1, a result 
not following from perturbation theory. It represents a transparent regime where 
three Bogoliubov states merge, leading to additional degeneracies of the topologically 
nontrivial ground state with odd number of Majorana fermions at the end of each wire. 
We also obtain the exact tunnel currents for a fixed parity of the eigenstates. The 
Josephson current shows the characteristic 47r periodicity expected for a topological 
Josephson effect. We discuss the additional features of the current associated with a 
closure of the energy gap between the energy levels. 
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1. Introduction 

In the Josephson effect [I] the phase coherent tunnehng across a junction between two 
systems A and B with broken ^7(1) symmetry imphes a dissipationless current oscillating 
with = — 0;,, where and 0b denote the phases of the superfiuid order parameter 
in A and B, respectively. The typical oscillating behavior in a conventional Josephson 
effect is given by / = Jq sin0, having periodicity 27r. Recently the Josephson effect has 
been considered in the framework of topological insulators |2]. It was shown that in 
a system where the tunnel junction is a topological insulator, a fractional Josephson 
current / oc ±sin(0/2) occurs as a consequence of the Z2 topological nature of the 
quantum spin Hall insulator. The occurrence of a contribution in the current featuring 
half of the phase difference has in this case a topological origin. In particular, the An 
periodicity arises due to boundary zero fermionic modes at the junction. Since these 
boundary modes have zero energy, they are not influenced by changes in the magnitude 
of the gap in the bulk superconducting state [21 E]. Interestingly, these zero energy 
boundary modes are found to be Majorana fermions, which in this context correspond 
to Bogoliubov quasi-particles having the reality property 7 = 7"'". In order to exhibit a 
concrete example of Bogoliubov quasi-particles fulfilling the Majorana reality condition, 
it is enough to recall the Bogoliubov transformation in superconductors [Ij, 

7lk = ^tkCk,a - ^kC^k,/3' 72k = ^ikC-k,a + ^^kcj,^^, (l) 

where ["UkP + I'^kP = 1- For conventional spin singlet superconductors we have a =t" 
and /3 =4, in which case the reality condition cannot be achieved. Thus, let us consider 
a = /3 or, more simply, spinless fermions. If Mk = and fk = e*^''/v^, we will 

have that at zero energy (i.e., /c = 0) 710 will be a Majorana fermion if ^0 = ^/2, 
while 720 will be a Majorana fermion if C,q = 0. In the specific model considered by 
Fu and Kane |21 E] , superconductivity was induced by proximity effect on the surface 
of a topological insulator, with surface excitations being described by a Dirac-like 
Hamiltonian in two dimensions. Within a mean-field treatment of the problem, the 
f/(l) symmetry breaking is introduced by coupling the anomalous fermion bilinears 
ip^ipl cind ipiip^ to an uniform superconducting gap A = Aqc^^. Thus, one can define 
spinless fermionic fields entering Eq. ([l) as Ck = {ipk,^ + e'^^'^ipKi) / For the case of 
two superconducting surfaces connected via a nanowire, one can assume the system to be 
one-dimensional and compute the corresponding fractional Josephson current [21 |5l |6] . 

As another instance of the theory considered in Refs. [2] and [5], it is interesting 
to notice that beyond mean-field theory its topological content in a nanowire setting is 
similar to the model of Goldstone and Wilczek in 1+1 dimensions [7] involving Dirac 
fermions coupled to two scalar fields (pi and (p2- Indeed, we can take (pi and as the 
real and imaginary parts of the superconducting order parameter, ip. In this way, the 
expectation value of the covariant fermionic current can be shown to be identical to 
the topological current [7], 

= if) = L,. A^*d.^-f.^*) _ (2) 
27r \(p\^ 
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Thus, if L is the length of the wire, the topological charge. 



Q = [ dx^K^ = — [ dx^ —— aictan ( — 
J 27r J_L/2 ox^ 



) 




— arctan ( — 



) 



L/2 



A(j) 
2^' 



(3) 



yields the charge of the topological excitation (soliton) in terms of the phase difference 
A0 between the boundaries only, and the order parameter in the bulk plays no role. 
Let us assume that each end of the wire is located in superconductors A and B, 
with corresponding phases (f)a and (pb, respectively. In this case we have ipi{—L/2) = 
Ao,aCos0a(-L/2), ip2{-L/2) = Ao,a siu 0a(-L/2) , (pi{L/2) = Ao,bCos0b(L/2), and 
{p2{L/2) = Ao,6sin(/)fc(L/2), so that 



which for both L/2) and 0b(L/2) in the interval (— 7r/2,7r/2) yields the principal 
value A0 = (f)h{L/2) — L/2). The function arctan(tan0) is double-valued at 
= {2k + l)vr/2. A; G Z and has principal value in the interval — 7r/2 < < 7i/2. 
Thus, as the phase at one of the ends of the wire vary from to 27r, the arctan(tan 0) 
makes an abrupt change of sign at tt/2. In the context of Eq. ^ the phases at the 
boundaries are constrained to half of the period of the order parameter, reflecting the 
fact that solitons (and Bogoliubov quasi-particles as well) are created in pairs, and the 
full period corresponds to two solitons. Therefore, we should take varying from — 7r/2 
to 7r/2, e. g., 0a(— L/2) = — 7r/2 and 0b(L/2) = n/2, so that A0 = vr, implying the 
fractional charge Q = 1/2. It is of course also possible to choose 0a (—L/2) = 7r/2 and 
0fc(L/2) = — 7r/2 and obtain Q = —1/2. This pair of solitons of charge 1/2 corresponds 
to the Majorana boundary modes in the wire. This fractional charge of a Majorana 
fermion can also be directly checked using the spinless variant of the Bogoliubov quasi- 
particles ([T]) at A; = 0. Indeed, using our previous discussion after Eq. ([T]), we see that 
for (q = tt/2 we have a particle number operator 7jo7io = 7io = 1/2, with a similar 
result holding for the Bogoliubov quasi-particle 720 when (q = 0. 

All the physical and mathematical properties discussed above are also present in 
a simple one-dimensional lattice model introduced by Kitaev some time ago [3]. This 
model has recently been used to motivate interesting realizations of physical systems 
where Majorana fermions may play a crucial role P, El |9l [ini ttH [ISl UHl [161 El UHl [19] • 
Within the framework of Kitaev's model, we will investigate two chains of spinless 
fermions having broken U{1) symmetry connected by a tunnel junction. We will first 
assume that the broken symmetry state is such that the superconducting gap |A| = w, 
where w is the hopping between nearest neighbor sites. This simplifying assumption has 
the advantage of allowing an exact derivation of the energy spectrum of the Josephson 
effect. The exact solution unveils interesting features at strong-coupling which may 
possibly be present in more realistic situations involving semiconducting wires having 



A0 = arctan[tan0;,(L/2)] — arctan[tan0a(— L/2)] 



(4) 
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strong spin-orbit coupling lying on the surface of an s-wave superconductor [9l[Tn]. One 
salient feature at strong Josephson coupling is the existence of a critical value of the 
tunnel coupling, Qc = \/2w, leading to level crossing between all the energy levels having 
different quantum numbers m, thus closing the gap between the levels and increasing 
the degeneracy. For such a value of the coupling the energy levels can be equivalently 
written as crossing levels of the form Ema{(p) = 2(Ta/2w cos(0/6 — 27r(m — l)/3), where 
m = —1,0, 1 and a = ±1, quite different from the leading order weak-coupling result 
featuring an energy ~ cos(0/2) [2], |9l [TO]. In order to check whether our results are 
generic, we also consider the numerical solution of the problem when |A| 7^ w and show 
that the same gap closing feature is present in the more general situation, although the 
structure of the spectrum is more complex. Despite the 0/6 factor in the argument 
of the cosine, the Josephson current still exhibits the An periodicity characteristic of 
the topological Josephson effect. In fact, we will show that the factor 0/6 actually 
corresponds to the principal value of a certain double- valued function ^(0) having 47r 
periodicity. At the same time, the current will be shown to exhibit some additional 
features which arise due to closure of the energy gap. 

The plan of the paper is as follows. In Section 2 we introduce the model and 
obtain the exact energy levels for the Josephson effect. Both situations corresponding 
to |A| = w and |A| ^ w are discussed. The former case will allow us to write down 
simple exact analytical expressions for the energy eigenvalues, while the latter case can 
be diagonalized exactly. In Section 3 we calculate the Josephson current from both 
equilibrium statistical mechanics and using an ensemble where the states have a fixed 
parity [20]. Section 4 concludes the paper. 

2. The Kitaev model based Josephson junction 

The Hamiltonian of the system consists of two wires A and B which are in contact via 
a tunneling Hamiltonian H = Ha + Hb + Ht [SI HD], where 

N-l 

Ha = ^[-w(ajaj-+i + a'^j^^aj) + Aaa^j+^a] + Alajaj+i], (5) 
i=i 

Hb = + bl,b,) + A,bl,b] + A*6,6,+i], (6) 

Hr = fi'(ajv&i + b\aN), (7) 

where A^ = lAle**^" and A^ = |A|e*'^'', with the magnitudes of the superconducting gaps 
assumed to be the same. If we perform the global gauge transformations aj — ?■ e^'^'^^'^aj 
and bj — )■ e'^'^''^%j^ the phases 0^ and 0?, are gauged away in the Hamiltonians Ha and 
Hbi while the tunnel Hamiltonian becomes, 

Ht = Tajv^i + T*b\aN, (8) 
where T = ge'^^'^, with = 0^, — 0^ being the phase difference across the junction. 
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We will first consider the case where |A| = w, where an exact analytic solution for 
the Josephson energy eigenstates is possible. 

It is convenient to write the Hamiltonians Ha and Hb in matrix form. Thus, we 



will write Ha = {1/2)iIj\MiIja, where ijj 



,t 



qn] and M is a 



symmetric 2N x 2N matrix having zero trace. The Hamiltonian Hb is written similarly 
with the Qj operators being replaced by the bj ones. Note that iPa^Pa = N. The matrix 
M has a doubly degenerated zero energy eigenvalue corresponding to two Majorana zero 
energy modes residing on each end of the chain, and (A^ — l)-fold degenerated energy 
eigenvalues ±2w. The Hamiltonian can be written in diagonal form by introducing the 
following fermionic operators, 
1 



71 



(ai + a[ 



In 



I'M 



V2' V2' 

which are the zero Majorana boundary modes, and the nonzero modes, 

1 ri, t t N 

-(-aai + aa[ + a2„+2 + a2n+2) 



(9) 



CN-(2n+l),a 



+ (1-5: 



V2n + 1 [2 

n 



m=l 



(10) 



1-5 



C7V-2(fc+l),cr 



N2 



1 , I I ^ 

-{ai — + a2k+i + «2fc+iy 



+ Y.^-ar^'am+2 



m=0 



:iii 



where a = ±1, n = 0, 1, ... , un, and k = 0, 1, ... , k^. Here = N/2 — 1 and 
kN = N/2 - 2 if is even, and = {N - l)/2 and kN = {N - 3)/2 if N is odd. The 



operators (|10l) and (|TT| correspond to two interpenetrating sublattices Li and L2. The 

(12) 



fermionic operators satisfy the local constraint. 



cr=±l 



which together with the constraint 

ihi = 7i = lilN = 1n = 1/2, 
for the zero modes yields 



7i7i + 7jv77V + XI 5Z S^-^J^- = 
i=i ^=±1 

Therefore, the Hamiltonian for a chain A with A^ sites has the diagonal form, 

N-l 

j = l a=±l 



(13) 



(14) 



(15) 



In view of the constraint (12) the above representation of the Hamiltonian Ha 



corresponds to localized spins 1/2 in an external magnetic field of magnitude 2w applied 
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along the 2;-direction. Note that we can use the constraint (12) in Eq. (15) to rewrite 
it in the form 

N-l 



Ha = 2w^l 



1/2), 



(16) 



obtained in Ref. PJ. 

The Hamiltonian Hb is diagonahzed in a similar fashion, except that we have to 
name the new fermionic operators differently, say dja^ for the nonzero modes, with Si and 
5n being the zero boundary modes. Since in the new basis the Hamiltonians Ha and 
Hb are written in terms of localized particles, it follows that the part of the spectrum 
corresponding to the Josephson energy levels can be obtained by solving a reduced 
Hamiltonian, basically one featuring two dimers {Ha and Hb for N = 2) connected by 
a tunnel junction. This procedure is basically the same as the one used by Alicea et al. 
[To] to solve the problem perturbatively. 

For the exact result we obtain a spectrum containing two zero energy eigenvalues 
corresponding to Majorana boundary modes, the energy levels ±2w, each with 
degeneracy 2(A^ — 2), and the phase-dependent Josephson energies obtained by solving 
the dimer/dimer junction, whose Hamiltonian is the same as H ioi N = 2, i.e., 

-ij^M^, (17) 



H 



where 



a\ al ai 02 b\ h\ ^1 ^2 



02 



a\ 
hi 
b2 
b\ 



M 






—w 





—w 














—w 





w 





r 














w 





w 














—w 





w 



















r* 











—w 





—w 














—w 





w 














-r 





w 





w 














—w 





w 






(19) 
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The above dimer/dimer problem is exactly diagonalizable by a Bogoliubov 
transformation. Before proceeding with this calculation, let us first diagonalize each 
dimer separately in order to gain some physical insight from the terms arising in the 
tunnel Hamiltonian. The dimer Hamiltonian for the first chain is diagonalized by the 
transformation, 

71 = ^(ai + 4), (20) 

72 = ^(4-«2), (21) 
= ^{ai- a\ + a2 + a\), (22) 

c+ = ^(-fli + aj + fla + 4)) (23) 

where 71 and 72 are Majorana boundary modes satisfying 7^ = 7I = 1/2, and the 
constraint 

cVc+ + c^_c_ = 1 (24) 

holds. Note that the equations above are just the general transformation given in Eqs. 
dOpOfnl ) for the special case N = 2. Thus, the Hamiltonian for the dimer A is in this 
new operator basis given by, 

Ha = w{c\c+-c^^c_) = 2w(^\c+-]^. (25) 

Note that the Majorana fermions do not appear in the Hamiltonian H^, since their 
energy eigenvalues vanish. We obtain a similar result for the Hamiltonian Hb-, with 
the new fermionic operators being called d+i and Here 5i and 82 are the 

corresponding Majorana boundary modes for the dimer Hamiltonian Hb- 

The Majorana fermion operators will appear explicitly only in the tunnel 
Hamiltonian connecting the states of the two dimers. This will make two of the four 
Majorana modes overlap, making in this way two zero modes disappear. This is better 
seen by collecting the Majorana fermions at the junction into a new (ordinary) fermionic 
operator defined by, 

/ = ^(72 + ^5i). (26) 

Thus, the Hamiltonian of the dimer/dimer system can be written in the form, 

H = Ho + H, (27) 

where 

Ho = w{clc+ - clc. + dld+ - did-) - g cos (^0 (^f^f - , (28) 
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and 



H = '-^[^icl + ct_)(/t - /) - ^(f + f^){d. - d+) 

+ c5_rf_ + c^_d- - c5_rf+ - c^_d+] + h.c. (29) 

The Hamiltonian H describes the hybridization between the Bogohubov quasi-particles 
from Hq and, in addition, tunnehng processes involving fused Majorana states across 
the junction. 

Let us first neglect the hybridization contribution and study the spectrum of the 
Hamiltonian Hq. The energy eigenstates are given by |n^j_, nl; ra^^, ral; rif), which are also 
eigenstates of the particle number operators A^^ = 4c±, Nf^ = did±, and Nf = f^f, 
since these operators commute with Hq. In view of the constraint (24), there are also 
corresponding constraints for the quantum numbers, i.e., 

nl + n'L = l, + 72^=1. (30) 

In Table 1 we show the eigenstates of Hq and the respective energy eigenvalues. We 
note that the eigenstates |1,0;0, 1;0) and |0,1;1,0;0) are twofold degenerate with 
energy eigenvalues -Eiooio = -E'oiioo = (^7/2) cos(0/2). A twofold degeneracy also 
occurs with the eigenstates |1, 0; 0, 1; 1) and |0, 1; 1, 0; 1), which have energy eigenvalues 
-f'looii = Equqi = —{g/2) cos(0/2). These degenerate energies vanish for = {2k + l)7i, 
with G Z. A closer look at Table 1 shows that we can attribute either a positive or 
negative sign to the coefficient of cos((/)/2), depending on whether or not the occupation 
number nj vanish. The coefficient of 2w is either +1, —1, or 0, depending on the 
configuration involving the occupation numbers n^. and n^. Thus, we see that the 
general expression for the energy eigenvalues of Ho can be written as, 

E^:i = 2mw +^ cos (^), (31) 



where 



m = -± + ~ -, a = (-1)"^ (32) 



with the understanding that the constraints (30) have to be satisfied, and Uf = 0,1. 
Thus, the quantum number m is physically the 2;-projection of the total pseudospin 
"magnetic" quantum number. The quantum number a = (— 1)"^ is the parity of the 
fused Majorana state. 

The exact energy spectrum is easily obtained from the diagonalization of the matrix 
(19), which leads to the secular equation, 

- + g^)E ± gw^ cos(0/2) = 0, (33) 

and the Hamiltonian of the effective dimer/dimer system describing the tunnel junction 
becomes, 

-^Junction ^ ^ -^mcr-^mcr; ("^^) 
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Table 1. Eigenstates |n^^, rtl; n^J., nl; ny) of Hq and their energy eigenvalues. The 
last slot in the eigenstates denote the presence or absence of a fused Majorana fcrmion. 



Eigenstate 


r ,11 f^vo"\r PI crpTi ^ri^ 1 1 1 p 

J-l/llCl pL y dclCli V cXl LLC 




1,0 


1,0 


0) 


-C/10100 — '^"J -|- 2 '-■Ub y^J 




0,1 


0,1 


0) 


-^01010 = -2W + 1 COS (1) 




1,0 


0,1 


0) 


-^10010 = f COS (1) 




0,1 


1,0 


0) 


-^01100 = f COS (1) 




1,0 


1,0 


1) 


^10101 = 2w - fcos (1) 




0,1 


0,1 


1) 


^01011 = -2w - fcos (1) 




1,0 


0,1 


1) 


^10011 = -f COS (f) 




0,1 


1,0 


1) 


^01101 = -f COS (f) 



where A^ma are the corresponding particle number operators. Therefore, we obtain along 
with a doubly degenerated zero energy eigenvalue, the six energy eigenvalues, 

nl/3 



-2i7r(m-l)/3 



a cos 



+ e 



2i7r(m-l)/3 



a cos 



+ iG(0) 



2 

nl/3 



- iG(0) 



where m 



-1, 0, 1 and (T = ±1 as before, and 



G(0) 



Aw"' 



Ag'^w'^ 

The following relation holds, 

along with the constraint, 

N^^ + Ar_„,_, = 1, 



cos^ 



(35) 

(36) 

(37) 
(38) 



as expected for Bogoliubov quasi-particles. This allows us to rewrite the effective 
dimer/dimer Hamiltonian (|34]) of the junction as, 

3 



H 



Junction 



i=i ^ ^ 



(39) 



where E\ = -Eo,-i) -E'2 = -E'l.+i, and i?3 = with a similar relabeling for the particle 

number operators. The form (39) will be useful in the calculation of the tunnel current 
later. 

The energy eigenvalues above have subtle analytic properties. In fact, in view of 
the property (37) and the oscillatory behavior of the energies, we can see that pairs of 
energy eigenvalues having the same m and opposite a necessarily cross at = (2/c + l)7r, 
k E Z. Thus, the in periodicity of the energies (35) is in this case equivalent to 27i- 
periodicity with double- valuedness at = {2k + l)7i, a fact already anticipated by Kitaev 
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[3]. Therefore, we can rewrite the energies (35) in terms of its principal values in the 
interval — vr < < vr as, 

27r(m-ir 



where 



g(0) 
_cos(0/2) 



^(0) 



arctan 



(40) 



(41) 



We emphasize that the cubic equations determining the exact Josephson energy 
levels always arise in the case where |A| = regardless of the size of the chains 
involved. For chains A and B having each lattice sites we obtain the characteristic 
polynomial for the corresponding 2N x 2N matrix M, 

P{E) = E\E^ - 4u^2)27V-4[^3 _ (4^2 ^ ^2)^ _ cos(0/2)] 

X [E'^ - (Aw^ + g^)E + gw^ cos(0/2)]. (42) 

Therefore, the Josephson energies are independent of the size of the wires involved, 
reflecting the topological nature of the Josephson effect in such a system. 
For g w we obtain up to second order in g/w, 



E ^ 



Aw ( 2TTm\ 

—pz sin 

V 3 y 

,2 

sin 



ag cos 



9' 



16^3 



27rm 



(5 



/ 27rm 
3 cos(/)), 



cos 



which yields. 



and 



ag 

2w + — cos 



ag 

—2w + — cos 



-ag cos 



+ w^i^ ~ 3 cos ( 



32w 



32w 



(5 — 3 cos ( 



(43) 

(44) 
(45) 

(46) 



agreeing with the perturbative result by Alicea et al. ^LO^. From the perturbative 
expansion it is seen how the degeneracy of the m = eigenstate of Ho is lifted when 



(f) 7^ {2k + l)7r, G Z. Note that while the first two terms in Eqs. (44) and (45) agree 
with the energy eigenvalues -Eioioo; -^'oioio; -E'loioi; and -Eoioii of Hq, Eq. (46) is twice 
the energies -Eiooio, -^oiioo, -Eiooii, and £^01100 (see Table 1). 



In Fig. [T] we plot the energies (35) for four different values of the ratio g/w, 
from weak- to strong-coupling. Note the crossing between the levels £'^(7(0) for a 
given m, with a gap between levels having different values of m. In panel (a) of 
Fig. [1] we have a typical weak-coupling situation, which may be well described by 
the approximate formulas, Eqs. (44), (45), and (46). However, as the strong-coupled 



regime is approached, the gap between the levels having different values of m starts to 
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Figure 1. Exact energy levels E„ia{<t>)/w for four different values of g/w: (a) 
5/w = \/2/2; (b) g/w = \/2/l-2; (c) g/w = ^2/1.09; (d) 5/w = \/2. The color 
scheme of the curves here refer to the following. The red and magenta curves are 
i?+i,CT(</>), the green and brown ones are £^_i_o-(</'), while the blue and the purple ones 
represent £^o,ct (</>)■ It is readily seen that curves Ema{<l>) for a given m cross, with a 
gap between the curves having different values of m. In Panel (d) the levels merge, 
being in this way equivalent to the energy curves given in Eq. (|52[). 



close [panels (c) and (d) in Fig. [l]. For g/w = a/2 [panel (d) in the Fig. [l] the gap closes 
completely, corresponding to a merging of three energy levels with quantum numbers 
m = —1, 0, 1. For g > y^w the gap opens again and it starts to grow for increasing g. 
For g < \/2w and for g > \/2w the levels Ema for a given m are doubly degenerated at 
= {2k + 1)11, A; G Z, while no degeneracies occur at = 2kn [See, e.g., panels (a), (b), 
and (c) at Fig. [l]. However, for g = \p2w double degeneracies also occur at = 2/c7r, 
with the difference that these degeneracies also represent points of non-analyticity of 
the energy eigenvalues as functions of 0. 

The level crossings at = vr are particularly important, especially the ones crossing 
zero energy. Indeed, the level crossings at zero energy are associated to additional 
Majorana zero energy modes which are a linear combination of fermionic operators 
living on both sides of the junction. Thus, these additional Majorana modes arising for 
= 7? are very different from the ones living at the outer boundaries (the left end of 
chain A and the right end of chain B), 71 and ^at, which are independent of 0. The zero 
energy Majorana modes at = tt are given by 
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72 



r2 + 4 



ir 



(aiv_i - ajv_ J + hi + h\ 



(4J 



where r = g/w. The Majorana fermion 71 is a superposition of the Majorana fermions 



In 



Oat 



and (^2 = {p2 + h\)/\/2, while 72 is given as a superposition of 



the Majorana fermions 7Ar-i = [i/ \/2){a}j^_-^ — a^^i) and 5i = (61 + b\)/\^. Thus, for 
(f) = n the total number of Majorana states (four, including the outer boundaries) is 
the same as the number of Majorana modes for two decoupled chains. Note that the 
decoupled chains limit r — )■ yields 71 = and 72 = 5i, as expected. These Majorana 
states at = TT are robust and continue to exist even at very large Josephson couplings, 
r — )• 00, in which case they become 71 = 62 and 72 = 7Ar_i. Since varying r corresponds 
to changes on the bulk properties of the system, due to the fact that |A| = w, the 
zero Majorana modes at = vr are topologically protected. Furthermore, their number 
= 7^ + 7I = 1, leading to a parity (—1)"^ = —1 at the tunnel junction when 



= vr. By recalling Eq. (38), we see that for m = and = vr the number operators 
^o,+i = 7i = 1/2 and A^o,-i = 72 = V^- Therefore, the charge is fractionalized in this 
case. 

The additional level crossings at g = \/2w arise because for this value of g the 
function G{(f)) in Eq. (36) becomes G(0) = | sin(0/2)|, which is non-analytic at = 2A;7r. 
The energy levels become, 

1/3 



-2m{m-l)/3 



a cos 



sm 



+ e 



27rj(m-l)/3 



cr cos 



+ i 



sm 



1/3" 



If we consider m = a = 1 and —2n < < 2tt, Eq. (49) simplifies to, 

Eii(0) = 2y2w;cos(0/6). 

The muultivaluedness of Ema{4>) leads to 

'0 27r^ 



En(0) = 2v^ 



w cos 



(49) 



(50) 



(51) 



6 3 

in the interval 2tt < (p < 6tt, which is the same functional form as -E'oi(0) in the interval 



< < 27r. By further analyzing the functional dependence of (49) on and its 



quantum numbers, we obtain that the spectrum for g = \/2w shown in panel (d) of Fig. 
[T]is indistinguishable from an energy spectrum of the form, 

27rfm - 11 



Emcr{4>) = 2cta/2wcos 



+ 



(52) 



6 3 

More precisely, the energy levels above are the principal values of the energies given in 



Eq. (49). This can be seen from Eq. (40), which for g = \/2w becomes 

27rfm - 11 



Emcr{4>) = 2a\/2w COS 



(53) 
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Figure 2. Double-valued function 9{(j)) = (1/3) arctan(| sin((/)/2)|/cos(0/2)) [Eq. 



(54)] 



where 



arctan 



|sm(0/2)| \ 
cos(0/2) J 



(54) 



lead to the curves shown in panel (d) of Fig. [T} Indeed, in Eq. (53) we can still see that 



is now a double- valued function of 0; see Fig. [2j Eq. (53) describes the nature of the 
special point g = \/2w more accurately than Eq. (52), although both sets of functions 



the 47r-periodicity of the spectrum via the double- valued function (54), while this is not 
any longer apparent when the functions (52) are used. Note that 6{(j)) is non-analytic 
at = 2kTc. From this we immediately see that for g = \/2w the Josephson current will 
exhibit jumps at = 2k7T. This expectation will be confirmed by explicit calculations 
in the next Section. 

For g = \/2w a tunneling event from the state having quantum numbers m = 1 
and cr = ±1 into one having m = and a = ±1 is allowed at = 2A;7r. In fact, 
we can see from panel (d) in Fig. [l] that there is for g = \plw an eigenstate with 
energy £"1 _i(0) = i?o,-i(0), thus allowing a transition to the eigenstate having quantum 
numbers m = and a = — 1, which crosses zero at = vr, corresponding to a Majorana 



zero mode there. Furthermore, if we use the energies (52), all curves cross zero for some 
value = (2k + l)7r. A similar behavior to Eq. (53) is obtained in a three-junction 
Josephson ring with a magnetic flux inside the triangular loop [23] • In this case m 



would correspond to the number of vortices. In the case of Eq. (53) each junction in the 



triangular loop would refer to a Z2 fractionalized Josephson effect, i.e., one featuring 
half of the phase difference. 



The result of Eq. (49 ) is quite interesting as it does not follow from the perturbation 



theory [10] and was not anticipated in the original work [3]. In order to see whether 
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this degeneracy found for g = \/2w is a general feature of the model we consider also 
the general case where |A| w. Although in this case an exact solution is still possible 
for each individual chain |3l El] , to solve exactly the Josephson system involving the 
two chains analytically is not an easy task. Observe that in this relevant case one has 
to take into account that chains must have free ends. As a result, the momenta will 
not belong to the first Brillouin zone, but will satisfy a transcendental equation [21] . 
Furthermore, zero boundary modes (Majorana modes) are only present for odd and, 
strictly speaking, the crossings of the energy levels at zero energy for = {2k + 1)tt do 
not any longer occur if A^ is not large enough. In Fig. |3] we show the results for the 
exact diagonalization of the Hamiltonian (l)-(3) for |A| = w/2 and 18 sites (A^ = 9 
for each of the chains). Fortunately, the value of A^ for which level crossing occurs in a 
similar way as in Fig. [T] is not too large and already for A^ = 9 results similar to the 
exact solution can be be found. In Fig. |3](a) we show the results of the diagonalization 
for g < gc = a/5/2. Observe that in contrast to the exact analytical results there 
is no actual crossing of levels at zero energy for = vr and A^ = 9, due to a small, 
almost negligible overlap of the Majorana fermions; see the inset in the Fig. [3]- (a). 
This is a finite-size effect and as A^ further increases, this small gap at = tt becomes 
exponentially small, so that the crossings at zero energy for = {2k + l)7r from Fig. [l] 
are approximately recovered. We have confirmed this behavior with system sizes up to 
A^ = 50 per chain. Another interesting result is that for increasing g to the critical value, 
which is here gc = y/5w/2, one finds additional crossings at = 2A;7r [Fig|3](b)], again in 
correspondence with the exact solution for the case where |A| = w. This confirms that 
the analytical results, obtained for |A| = w, are quite generic and do not depend on the 
choice of the parameters. 



3. Topological Josephson current 

The standard way for calculating the Josephson current is given by the following formula 
(in unities of 2e/h), 

where -F(0) is the free energy of the junction as a function of the phase difference 0. 
In the above equation all states are taken into account, since standard equilibrium 
thermodynamics is used, i.e., the partition function is simply given hj Z = 
Tr[e~^^'^~^^-']. However, in the case studied here the parity of the states involved 
plays an important role. This point has been extensively discussed in the literature 
[21 El [m [131 E] • Superconducting tunneling in systems with fixed parity [1] grew in 
importance in the 1990s in view of experiments performed with small tunnel junctions 
[20] . In the case of topological superconductors the role of fermionic parity is even more 
crucial, as the parity of states is associated to the presence or absence of Majorana 
boundary modes [3l [H] . 




Figure 3. Energy levels for the tunnel junction involving two Kitaev chains, each 
with = 9 (making a total system size of 18 sites), and |A| = w/2, for (a) g < gc 
and (b) g — gc = y/5'w/2. The inset in (a) zooms the region around = tt for g < ge- 
mote that there is actually no exact crossing at = tt and at zero energy, in contrast 
to |A| = w, which indicates a small overlap of Majorana fermions. 



The correct tunnel current is calculated by projecting out fixed parity states in the 
partition function. In order to do so, we first note that only three states are physically 
distinct. Indeed, as usual in any Bogoliubov treatment of superfiuid systems, redundant 
states are introduced and quasi-particles occur in pairs and the energy spectrum includes 
in this way the energies of particle and hole states [Ij. In our case, this fact is expressed 
in relation (37). Therefore, we can compute the partition function for fixed parity using 
just three states out of the six ones determined by Eq. (35), corresponding to the energy 
levels El = Eq _i, E2 = and E^ = -Ei _i. In other words, we simply have to use 

the effective Hamiltonian (39). If = A^^i + A^2 + ^3, the projection operators for even 
and odd parity states are given by [24j , 

^e/o=^[l±(-in. (56) 

The partition functions for fixed even and odd parities are 

Ze/o = Tr[Pe/oe-^''^--«°"], (57) 

respectively. Thus, in Ze only even parity states are retained, with odd parity states 
being suppressed. The opposite is true for Zo, where even parity states are suppressed. 
For the computation of the Josephson current we have to use the free energy 



m = b^(^]- (5J 



/3 \Zo 

In Fig. |4] we show plots of the Josephson current in the case of fixed parity and low 
temperature. The plots shown in Fig. |4] are essentially zero temperature ones; there 
is no appreciable difference between the current profiles shown and the ground state 
result. 
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Figure 4. Plots showing the Josephson current for fixed parity and low temperature 
(T = w/10). In Panel (a) we show the Josephson current for two different values of g, 
namely, g/w = \/2/2 (solid line) and g/w = a/2/4 (dashed line). Panel (b) shows the 
Josephson current for the critical value of the Josephson coupling, g/w = \/2. 



Nevertheless the closure of the gap between the different energy levels aX g/w = \/2 
remains manifest in the calculations of the current at fixed parity. This is clearly 
seen from FigQb) where discontinuous jumps at = 2'Kk occur due to additional 
degeneracies in the system. Mathematically the origin of the discontinuities comes 
from the fact that the function in Eq. (36) becomes G(0) = |sin(0/2)| for 

g = \plw^ which is non-analytic at the points 0^ = 2A;7r, making dG/d(^ discontinuous 
at these points. We can interpret these jumps physically in a way similar to the one 
of a so called SNS junction, where the energy of an Andreev bound state is given by 
E = Aa/1 — sin^(</)/2) (we have considered for simplicity only one bound state) [25] . 
In the case of a fully transparent normal interface the transmission coefficient becomes 
T = 1, and E = A|cos(0/2)|. Thus, inspired by this result, we can interpret the 
quantity, 

^ ~ ~/ ^^72' ^^^^ 

as a transmission coefficient. For g = \/2w the junction becomes fully transparent. 

We also note in passing that the exact results at weak-coupling, although being 
qualitatively similar to the second-order perturbative result [10] , still differ with respect 
to the the maximum value of the current, as shown in Fig. [5j At the same time, for 
the critical g/w = \/2 the perturbative results look completely different as they miss an 
additional feature associated with the closure of the energy gap. 

It is well known that superconductivity cannot occur in one-dimensional systems 
at finite temperature. However, in the Kitaev model the U{1) symmetry is explicitly 
broken by proximity effect. We are not dealing with spontaneous symmetry breaking 
here, which is actually a property of the higher dimensional substrate over which 
the wire is placed. Indeed, the gap in the wire is held fixed, a situation that may 
be approximately achieved using a proximity effect with a superconductor at higher 
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Figure 5. Comparison between the exact result (solid line) and the second-order 
perturbative result [TU] (dashed line) for g = a/2/2 and T jw = 0.01. In this example 
the perturbative result overestimate the maximum current by about 11 %. 






Figure 6. Plots showing the Josephson current for fixed parity and high temperatures. 
In Panel (a) we show the Josephson current for gjw = \/2/2 and two different 
temperatures, namely, Tjw = 2 (solid line) and Tjw — 1 (dashed line). Panel (b) 
shows the Josephson current for the critical value of the Josephson coupling, gjw — \/2 
and Tjw = 2. 

dimensionality. The wire can be assumed to be made of a semiconducting material 
with strong spin-orbit coupling, which eventually becomes an one-dimensional p-wave 
superconductor via proximity effect with the surface of an s-wave superconductor [3l HD] • 
In Fig. |6] we show plots of the Josephson current at finite temperatures. Not 
surprisingly, the amplitude of the Josephson current decreases with the temperature [see 
Panel (a)], while the An periodicity remains intact. Observe also that the Josephson 
current for g/w = \/2 again shows a discontinuous behavior which points towards 
the special character of the spectrum at this value of g/w. Observe that here the 
discontinuous jumps are 27i modulated just like in the low temperature case, but the 47r 
periodicity remains intact as it should. 
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4. Conclusion 

We have analyzed the Josephson effect for two lattice quantum wires featuring fused 
Majorana zero energy boundary modes at the tunnel junction. In the weak-coupling 
regime the exact solution reproduces the perturbative result [TU] for the energy 
containing a contribution ~ ±cos(0/2) relative to the tunneling of paired Majorana 
fermions. As the tunnel amplitude g grows relative to the hopping amplitude w, the gap 
between the energy levels gradually diminishes until it closes completely at the critical 
value (7c, whose magnitude depends on the ratio of |A|/iy. For |A|/w = 1 and gc = \/2w, 
the Josephson energies can be cast in the form given by Eq. (49), which is very different 
from the result obtained at weak-coupling. Although this regime which occurs at g > w 
is rather exotic from the point of view of its experimental realization, it is still interesting 
to see that Kitaev's model is richer than it was originally anticipated. In addition, the 
experimental setup for the realization of the Majorana fermions in quantum wires is 
still under discussion; see for example Ref. [B]. Thus, it would be interesting to analyze 
whether an experimental setup allowing for large values of the coupling g between the 
chains can be engineered. Furthermore, such a system can in principle be engineered 
using ultracold fermions in an one-dimensional optical lattice [26] . In this case a way to 
achieve the strongly coupled limit would be to consider the distance between the chains 
as being smaller than the lattice spacing in the chains. On the other side, it should 
be noted that the closure of the gap at gc may be spoiled in several ways in realistic 
systems, since additional Andreev states occuring at larger Josephson couplings may 
lead to a small gap. Even in the context of the simple fine-tuned model we solved here 
something different may occur if the quantum character of the phase difference, which 
is important for sufficiently small systems, is taken into account. By this we mean to 
consider the phase of the order parameter as a quantum operator conjugate to particle 
number. In spite of the difficulties involving the definition of a hermitian phase operator 
[27j . a semi-classical treatment is possible in Josephson nanosystems and the gap closing 
effect we have found may disappear due to such a contribution [23] • However, in a setup 
where superconductivity is induced via proximity effect, the simple description in terms 
of the Kitaev model may apply, and the new aspects we have discussed here is likely to 
be remain robust . Another interesting question is whether the additional degeneracies 
remain intact if interaction effects are included in the Kitaev Hamiltonian [12], or if 
disorder is present in the system. As far as the topological stability of the Majorana 
fermions at the boundaries is concerned, recent work indicates that they are stable 
against disorder [2S]. This does not necessarily mean that our additional degeneracies 
and 1/6 fractionalization ai g = gc survives disorder effects. This aspect of the problem 
needs further investigation, being beyond the scope of the present study. 

One further interesting aspect of the problem we have discussed concerns parity 
effects on the Josephson current. The standard equilibrium calculation where the parity 
is allowed to change would completely fail to account for the 47r periodicity of the 
topological Josephson effect. In the standard equilibrium calculation a sudden change 
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of sign would occur in the Josephson current when cj) = {2k + 1)71, which in turn would 
spoil the Air periodicity in the current, despite the An periodicity of the energies. On 
the other hand, when the calculation is done using a fixed parity ensemble |20], the 
discontinuity jumps associated with the transition between the states with different 
parity disappear and the 47r periodicity is restored. Nevertheless the closure of the gap 
between the different energy levels at g/w = a/2 shows up via discontinuities in the 
current as a function of the phase, but these discontinuities are not of the same type as 
the ones appearing in the calculations where the parity is allowed to change. 
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